Transition Temperature of the Homogeneous, Weakly Interacting Bose gas 
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We present a Monte Carlo calculation for up to N ~ 20 000 bosons in 3 D to determine the shift 
of the transition temperature due to small interactions a. We generate independent configurations 
of the ideal gas. At finite N, the superfluid density changes by a certain correlation function in the 
limit a — > 0; the N — + oo limit is taken afterwards. We argue that our result is independent of the 
order of limits. Detailed knowledge of the non-interacting system for finite JV allows us to avoid 
finite-size scaling assumptions. 
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Feynman Q| has provided us with a classic formula 
for the partition function of the canonical noninteracting 
Bose gas. It represents a "path-integral without paths" , 
as they have been integrated out. What remains is the 
memory of the cyclic structure of the permutations that 
were needed to satisfy bosonic statistics: 



Zn= ^({ m fc}); withP({m fe }) = Yl 



rrik 

Pk 



{m k } 



k=l 



mk\k mk 



(1) 



The partitions {m-fc} in eq. (|T|) decompose permutations 
of the N particles into exchange cycles (m., cycles of 
length i for all 1 < i < N with J2k ^ m k — N). pk 
is a system-dependent weight for cycles of length k. 

In this paper we present an explicit Monte Carlo cal- 
culation for up to ~ 20 000 bosons in three dimensions, 
starting from eq. ([!]) . The calculation allows us to deter- 
mine unambiguously the shift in the transition temper- 
ature T c for weakly interacting bosons in the thermody- 
namic limit for an infinitesimal s-wave scattering length 
a. This fundamental question has lead to quite a num- 
ber of different and contradictory theoretical as well as 
computational answers (cf, e.g., [0-0]). 

We will first use Eq. (|l|) and its generalizations to de- 
termine very detailed properties of the finite- N canoni- 
cal Bose gas in a box with periodic boundary conditions. 
We then point out that all information on the shift of 
T c for weakly interacting gases is already contained in 
the noninteracting system. In the linear response regime 
(infinitesimal interaction) , it is a certain correlation func- 
tion of the noninteracting system which determines the 
shift in T c . This correlation is much too complicated to 
be calculated directly, but we can sample it, even for 
very large N. To do so, we generate independent bosonic 
configurations in the canonical ensemble. We have found 
a solution (based on Feynman's formula Eq. (|l|)) which 



avoids Markov chain Monte Carlo methods. In our two- 
step procedure, a partition {rrik} is generated with the 
correct probability V({rrik}). Then, a random boson con- 
figuration is constructed for the given partition. 

We stress that all our calculations are done very close 
to T c , so that the correlation length £ of any macroscopic 
sample is much larger than the actual system size L of 
the simulation. This condition L -C £ allows us to invoke 
the standard finite-size scaling hypothesis || , but also to 
take the N — > oo limit after the limit a — ► 0. 

A key concept in the path integral representation of 
bosons is that of a winding number. Consider first the 
density matrix p(r, r', (3) of a single particle [r = (x, y, z)] 
at inverse temperature f3 = 1/T. In a three-dimensional 
cubic box of length L with periodic boundary conditions, 
p(r, r', (3) = p(x, x',f3)x p(y, y', (3) x p(z, z', (3) with, e.g., 



p(x,x',/3)= 



exp[-(x- {x 1 + Lw x )) 2 /2(3] 



(2) 



w x — — oo 



In Eq. (0), x and x' are to be taken within the periodic 
box (0 < x, x' < L). 
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FIG. 1. A one-dimensional periodic simulation box with 
three bosons. Particles move in imaginary time < r < j3 
and in periodic space < x < L. We use non-periodic coor- 
dinates. 
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It is more convenient to adopt non-periodic coordinates 
(—00 < x,x' < 00), as we will do from here on. In Fig. 
1, the path drawn with a thick line can equivalently be 
tagged by (xi,x[) or by (xi,x[). This notation allows 
one to keep track of the topology of paths without in- 
troducing intermediate time steps r, even for very small 
systems. With this convention, the winding number of a 
configuration, W = (W x , W y , W z ), is defined as 



w = 5>< - n )/L. 



(3) 



The winding number W in Eq. (^) is the sum of the (in- 
teger) winding numbers for each of the statistically un- 
corrected cycles which comprise the configuration. The 
complete statistical weight pk of a cycle of length k [cf. 
Eq- (0)] is given by the sum of the weights pk, w for all 
winding numbers w: 
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Pollock and Ceperley M have obtained the result 

< W 2 > L 2 
Ps/P= 3f3N ' (5) 

which connects the system's superfluid density p s /p to 
the winding number in a rigorous fashion. 

It is possible to determine the mean square winding 
number < W 2 > from Eq. (|j). We first compute < W 2 > 
for a given partition {m^} 



< w 2 >{„ lfc }= < w 2 > k 



(6) 



Here, < w 2 >k is the mean with respect to cycles of 
length k, < w 2 > fe = 3E„ Pk,wW 2 ][J2 w Pk,w? / Pk- T nis 
yields, by summation over partitions 



<ff2 >= El <w2 >fc Z N-k/Z N . 



(7) 



We have also determined the probability distribution of 
W x . 

An analogous calculation formally replaces < w 2 >k~^ 
k in Eq. (0), which becomes ^ fc k mk — N [cf. Eq. ([!])]. 
Equation (fa) is transformed into 



Zn = X Pk Z N-k/N. 



(8) 



Equation <ffl) allows the recursive calculation of the par- 
tition function Z]$ if Z\, . . . , Zisr-% are known 0. 
The same relation Eq. @ allows us to identify 



k < m k >= pkZ N - k /Z 



N 



(9) 



as the mean number of particles in a cycle of length k. 

From a different point of view, the quantity 
~Yj k kmke~P kei I pk determines the occupation number of 



single-particle energy levels q for a given partition {mk}- 
This allows us to compute the average number < iVj > 
of particles occupying state e, in the bosonic system: 



<Ni >= 



fe=l 



J N 



(10) 



Eq. ( ]To| ) is of crucial importance: We find that Nq/N, 
the condensate fraction, is different from the superfluid 
fraction, as determined from Eqs (|5|) and (|7]) in a finite 
non-interacting system. 

The term {} in Eq. ([l(]) can be regarded as the proba- 
bility P(n,i > k) of having at least k particles in state e^. 
Taking the sum over all states i, with the use of Eq. (||), 
we can connect cycle statistics with the usual occupation 
number representation: 



k < m k >= ^ p ( n i > k )- 



(11) 



This curious result, which is of practical use in inhomo- 
geneous systems [||, tells us that the discrete derivative 
of the mean cycle numbers with respect to their length is 
given by the probability of having k particles in the same 
single-particle energy level. 

Rescaled superfluid densities N^ 3 p s /p (from Eqs (]?]) 
and @) are plotted in Fig. 2a for Ni = 37, N 2 = 
296, N 3 = 2368, N 4 = 18 944 as a function of the rescaled 
temperature t = (T — T£°)/T£°, where T c °° is the critical 
temperature for N — ► 00. 
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FIG. 2. Rescaled superfluid density of an ideal Bose gas. 
The curves for different N with Ni+i — 8Ni intersect approx- 
imately at T c °° , as shown in (a) . The close-up view (b) reveals 
important differences. We determine the shift of the intersec- 
tion points as a function of the interaction (light arrows) . The 
dark arrow shows schematically the extrapolated shift in the 
thermodynamic limit. 

A finite-size scaling ansatz, which was used in previ- 
ous Monte Carlo work on the problem || , assumes that 
the curves of N^ 3 p s /p for a weakly interacting Bose gas 
should intersect at the transition temperature, as they do 
approximately. However, the small-scale Fig, 2b clearly 
shows the importance of corrections to scaling (cf. jt|) 
already for the noninteracting gas. By continuity, the 
corrections to scaling for the weakly interacting Bose gas 
must be important, especially if the temperature shift 
due to interactions becomes small. 
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Our strategy greatly benefits from the solution Eq. (^) 
for the ideal gas. We compute the intersection point 

1 /3 

(N i p s /p,t) for two finite systems with Ni and N2 = 
8iVi particles and determine how this point is shifted 
under the influence of interactions (cf. Fig. 2b). Our 
arbitrary but fixed ratio N2 /N\ = 8 facilitates the direct 
extrapolation in N± to N\ — ► 00. 

To generate a random partition, we interpret the term 
PkZN-k/Z]y in Eq. (|J) as the probability to split off a 
cycle of length k from a configuration of N bosons, and 
to be left with a system of N — k bosons. We can pick k 
with probability ~ pkZ[j-k with a simple "tower of prob- 
abilities" strategy [jl(| . Recursively, we can thus generate 
an independent random partition {m*,} with great speed. 
The recursion stops as soon as we have split off a cycle 
of length j from a system with j particles. 

To go from a random partition to a random configura- 
tion, we may treat each cycle separately. For a cycle of 
length fc, we select a winding number w x with probabil- 
ity Pk,w x [cf- Eq. (0)], and analogously for w y and w z . 
Towers of probabilities are again used. The cycle starts 
at a random position r = (x,y,z) with < x,y,z < L, 
and ends at r' = (x + w x L,y + w y L,z + w z L). Inter- 
mediate points are filled in with the appropriate Levy 
construction We have tested our algorithm success- 
fully against the known results (cf. Fig. 2). 

We thus generate independent free boson path-integral 
configurations by a method very different from what is 
usually done in path-integral (Markov-chain) Quantum 
Monte Carlo calculations, but with an equivalent out- 
come: Any appropriate operator is sampled with the 
probability 
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Here, indicates the summation over all permutations 
P, and r[ is the position of the particle P(i). In the pres- 
ence of interactions, the statistical weight of each con- 
figuration is no longer given by the product of the one- 
particle density matrices no = [Yii p( r iT r 'i)]- To lowest 
order in the interaction, the density matrix is exclusively 
modified by s-wave scattering. Likewise, only binary col- 
lisions need to be kept. This means that the correct 
statistical weight is given by 

TT a {r 1 , . . .,r N ;r[, ...,r' N ) = tt JJffy (n, r jt r-, r'j). (13) 

The contribution of collisions is to lowest order in a 

Y[ 9ij = l-aJ2cij (14) 
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Here, nj = ri — rj and 7y is the angle between and r\y 
Equation (|l^) corresponds to the popular Path-Integral 
Monte Carlo "action" with the following important mod- 
ifications: i) no interior time-slices are needed, ii) the 
interaction may be treated on the s-wave level, and Hi) 
the interaction may be expanded in a. For a consistent 
evaluation of the interaction with periodic boundary con- 
ditions, as schematically represented in Fig. 1, it is best 
to sum over all pairs i < j shown, with the condition that 
Ti be in the original simulation box (shaded in gray). As 
indicated by the small circles in Fig. 1, may have 
important contributions stemming from more than one 
representative of the path (rj ,r'j), especially for small 
systems. Of course, a cutoff procedure can be installed. 

We now find for the mean-square winding number in 
the interacting system 



<w2>a =<(l-aC)W>> 



< (1-aC) > 



(15) 



where we put C — Yli<j c v- Expanding in a, this yields 
< W 2 > a — < W 2 > = -a < (AW 2 )(AC) > , (16) 

where AO = O- < O > @- 

For a finite system of N bosons, the shift in the super- 
fluid density Sp s — p s {a) ~ p s (0) can thus be proven to 
be linear in a 



Sp s 
P 



X 



(3N 1 / 3 



(17) 



with X N =< (AW 2 )(AC) > /(3p). 

To determine quantitatively the shift of the intersec- 
tion points in Fig. 2, we expand the ideal gas superfluid 
density around the intersection temperature T s of two 
systems with N and 87V bosons at the same density: 

pJpi^N 1 / 3 = p s /p{T s )N x ^ + a N x[T- T s ). (18) 

In this formula, the linear expansion coefficients can be 
computed. With interactions, only p s j ' piT^N 1 / 3 is mod- 
ified to linear order in a. a(N) remains unchanged, as 
we restrict the expansion to \T — T s \/T s ~ ap 1 / 3 . We 
find the new intersection point of the two systems to be 
shifted in temperature as 



T s (0) 



T s (a)-T s (0) X ; 



8N 



X 



Ts(0) 
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(19) 



Iry-^expOr^ll^Kl + cos^)^]. 



We have also computed the shift in p s /p, but found out 
only that it must be extremely small. We are unaware of 
any fundamental reason for a vanishing shift in this quan- 
tity. In Fig. 3, we plot the shift AT S for different system 
sizes ranging from (Ni,N 2 ) = (37,296) to (2368,18944) 

— 1/2 

vs N x . We have not attempted a thorough analysis 
of the finite-size effects, which already appear negligible 
for our largest systems. 
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FIG. 3. 
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Shift of the intersection temperature 
ap 1/z ] as a function of N~ 1/2 (N 2 = 8AT X ). The 
are {Ni,N 2 ) = (37, 296), (125, 1000), (296, 2368), 
, and (2368, 18 944). 
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We conclude that the transition temperature of the 
weakly interacting Bose gas increases linearly in the scat- 
tering length a by an amount of 



(2.3 ± 0.25)ap 



1/3 



(20) 



Our result Eq. ( p0| ) is almost an order of magnitude 
larger than what was found in a previous Monte Carlo 
calculation ||. However, this calculation was restricted 
to very small particle numbers and it used a problematic 
finite-size scaling ansatz, as pointed out. The agreement 
of Eq. ( pfj| ) with the renormalization group calculation j^] 
seems to be quite good. 

It is very interesting to understand whether the result 
Eq. ( p0| ) directly applies to the current Bose-Einstein con- 
densation experiments (cf, e.g. jl2],[l3]] ) . In earlier papers 
H[l4), we have pointed out the particularities of these 
finite systems in external potentials (cf. |l5) for a gen- 
eral overview). Notwithstanding the differences between 
the two systems, a relevant parameter is for both cases 
ap 1 / 3 , where the maximum density (at the center of the 
trap) at the transition point must be taken in the inho- 
mogeneous case. The experimental value is of the order 
ap 1 ' 3 ~ 0.02. 

Within our method, we can also study finite values 
of the interaction, even though we no longer compute a 
correlation function, and also have to introduce interior 
time slices. Contributions beyond s-wave scattering need 
to be monitored, as we have in For N = 125 bosons 
we have found agreement with the linear reponse formula 
Eq. © up to ap 1 / 3 < 0.005, but a 15% decrease for the 
full treatment for ap 1 / 3 — 0.023. A detailed investigation 
of this question goes beyond the scope of this paper. 

In conclusion, it is worth noting that we have encoun- 
tered none of the difficulties which usually haunt boson 
calculations: We work in the canonical ensemble; there- 
fore, the fluctuation anomaly of the grand-canonical Bose 
gas plays no role. The density remains automatically con- 
stant as a function of a so that an expansion in ap 1 ! 3 is 
well-defined. At finite N, we can furthermore prove that 



the shift in p s /p is linear in the interaction parameter. 
We have also consistently approached the weakly inter- 
acting system from the vantage point of the ideal gas. 
This allows us to obtain the crucial information on ex- 
actly where to do our simulation (cf. Fig. 2). Finally, 
our extremely powerful sampling algorithm has allowed 
us to partially dispel the curse of Monte Carlo simula- 
tions: limitations to small system sizes. 
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